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Abstract. During the epoch of reionization, Ly-a photons emitted by the first stars can couple the neutral hydrogen 
spin temperature to the kinetic gas temperature, providing the opportunity to observe the gas in emission or 
absorption in the 21-cm line. Given the bright foregrounds, it is of prime importance to determine precisely the 
fluctuations signature of the signal, to be able to extract it by its correlation power. LICORICE is a Monte-Carlo 
radiative transfer code, coupled to the dynamics via an adaptative Tree-SPH code. We present here the Ly-a part 
of the implementation, and validate it through three classical tests. Contrary to previous works, we do not assume 
that Pa, the number of scatterings of Ly-a photons per atom per second, is proportional to the Ly-a background 
flux, but take into account the scatterings in the Ly-a line wings. The latter have the effect to steepen the radial 
profile of Pa around each source, and re-inforce the contrast of the fluctuations. In the particular geometry 
of cosmic filaments of baryonic matter, Ly — a photons are scattered out of the filament, and the large scale 
structure of Pa is significantly anisotropic. This could have strong implications for the possible detection of the 
21-cm signal. 
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1. Introduction 

The Epoch of Reionization (EoR) extends from the time 
when the first sources form as a result of the nonlinear 
growth of primordial density fluctuations in a fully neu- 
tral universe (z ~ 20), to the moment when the inter- 
galactic medium is fully reionized under the effect of UV 
radiations emitted by the sources (z ~ 6). This simple 
picture, however, hides a number of uncertainties. At this 
time, there are only two observational constraints on the 
EoR. The first comes from the Gunn-Peterson effect in 
the absorption spectrum of high redshift quasars. Indeed 
the transmitted flux drops sharply as a small neutral frac- 
tion appears toward high rcdshifts: it implies an ionization 
fraction xm < 10~ 4 at z < 5.5 (Fan et al. 2006). The sec- 
ond constraint is set by the Thomson scattering of CMB 
photons by the free electrons produced during the EoR. 
The corresponding optical depth, r = 0.09 ±0.03 (Spergel 
et al. 2007), implies that the intergalactic medium is al- 
ready significantly reionized at z — 11. In the next decade, 
we will learn a lot more by direct observation of the red- 
shifted 21-cm emission from the neutral IGM during the 
EoR. Such instruments as LOFAR, PAST or MWA will 
be able to probe the statistical properties of the 21-cm 
signal, while SKA will be able to make a full tomography 
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of the IGM up to z ~ 11. See Cariili et al. (2004) and 
Carilli (2006) for detailed prospects. 

In the last decade, a lot of work has been done, the- 
oretical and numerical, to predict the properties of the 
21-cm emission and better optimize the design of the fu- 
ture instruments. Madau et al. (1997) and Tozzi et al. 
(2000) present the first theoretical models of 21-cm emis- 
sion. The signal can be seen either in emission or in absorp- 
tion against the CMB, depending on the spin temperature 
of the neutral hydrogen. Interaction with CMB photons 
couples the spin temperature to the CMB temperature in 
less than 10 5 years during the EoR, which would make the 
signal undetectable. Fortunately two other processes tend 
to couple the spin temperature to the gas kinetic temper- 
ature instead. The first is the excitation of the hyperfme 
transition through collisions with electrons or other hy- 
drogen atoms (see Kuhlen et al. 2006 for numerical sim- 
ulations), which, however, is efficient only in overdense 
regions (baryonic Sp/p > 5/[(l + z)/20] 2 ). The second 
is the pumping of the 21-cm transition by Ly-a photons 
(Wouthuysen-Field effect, Wouthuysen 1952, Field 1959) 
which requires a threshold value for the local Ly-a flux to 
be effective. Consequently, the value of the kinetic tem- 
perature of the gas relative to the CMB temperature is 
crucial for determining the 21-cm emission brightness tem- 
perature. During the EoR, the gas is cooling down due to 
the expansion of the universe faster than the CMB, but is 
heated by hydrodynamical shocks from structure forma- 
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tion, X-ray from pop-Ill star or quasars and, to a much 
lesser extend, by Ly-a photons (see e.g., Furlanetto et al. 
2006a). Simple analytical models are now available which 
take into account various source types and formation his- 
tory (Furlanetto 2006b). They usually predict that the 
signal can be seen in absorption early on, then in emis- 
sion later on. The prediction for the typical amplitude of 
the differential brightness temperature is a few 10 mK. 
Analytical models cannot, however, take into account the 
full complexity of the 3D inhomogeneous IGM: numerical 
simulations including the dynamics of structure formation 
and, usually as a post-treatment, a full 3D radiative trans- 
fer are required. 

Dynamical simulations of structure formation have a 
long history, but cosmological radiative transfer simula- 
tions are a more recent field of investigation. Simulations 
showed that the ionized bubbles around the first sources 
are not spherical: indeed the ionization fronts propagate 
fast in the void and more slowly in the high density fila- 
ments (see Abel et al. 1999 for the first simulation). The 
geometry of reionization is now studied in large simulation 
boxes (~ 100 Mpc) to get a statistical sample of ionized 
bubbles of various sizes (Iliev et al. 2006a, McQuinn et 
al. 2007). However, with such large boxes, even with very 
high resolution simulations, small scale density structures, 
also called minihaloes, are not resolved and are often taken 
into account through a simple clumping factor. The global 
effect of minihaloes is to slow down reionization by con- 
suming photons during their photoevaporation. But a sim- 
ple, uniform, clumping factor may be insufficient to model 
the role of minihaloes, and more detailed studies have been 
performed (Ciardi et al. 2006, Iliev et al. 2005, and Shapiro 
et al. 2006 for the impact on the 21-cm signal). Using the 
density field from the dynamical simulations and the ion- 
ization fraction and gas temperature from the radiative 
transfer simulations, it is possible to produce 21-cm emis- 
sion maps. 

A number of authors provide predictions for different 
aspects of the 21-cm signal: the emission map at a given 
redshift, the average signal as a function of redshift or 
the signal power as a function of the angular scale. These 
quantities were first predicted from rather small simula- 
tion boxes (~ 10 Mpc) limiting the angular information to 
9 < 10 arcmin (Ciardi & Madau 2003, Gncdin & Shaver 
2004, Furlanetto et al. 2004, Valdes et al. 2006). Recently, 
predictions from larger simulation boxes (~ 100 Mpc) be- 
came available (Mellema et al. 2006, Iliev et al. 2007). It is 
a fact that the predictions, in particular the duration and 
intensity of the absorption phase, depend crucially on the 
modeling of the sources. But other factors have the poten- 
tial to alter these predictions. Indeed, all predictions from 
simulations, at this time, assume a uniform efficiency for 
the Wouthuysen-Field effect. However, Barkana & Loeb 
(2005) or Furlanetto et al. (2006a) recognized that fluctu- 
ations in the local Ly-a flux can produce additional fluctu- 
ations in the brightness temperature. An accurate quan- 
titative modelization of these fluctuations is important. 
Indeed, the 21-cm signal will be difficult to detect due to 



brighter foregrounds: the unique signature of its brightness 
fluctuations will make the extraction possible. Barkana & 
Loeb used simplified analytic models (neglecting radia- 
tive transfer effects on the local Ly-a flux) to compute 
this contribution. Our goal in this paper is to investigate 
how computing the full radiative transfer in the Ly-a line 
in a cosmological inhomogeneous medium can modify the 
picture presented by these authors. The 21cm signal will 
be difficult to detect, due to brighter foregrounds, and the 
unique character of its brightness fluctuations is a signa- 
ture that will make its extraction possible. It is therefore 
crucial to precise in more details those fluctuations. 

Although Monte-Carlo simulations of Ly-a transfer 
have a long history, starting with Avery & House (1968), 
only recently has the computing power become sufficient 
to tackle the case of a 3D inhomogeneous medium with 
kinematics, without restrictions on the optical thickness 
regime (Ahn et al. 2001). Several authors have now devel- 
oped similar codes to simulate the Ly-a emission from 
high redshift galaxies (Zheng & Miralda-Escude 2002, 
Cantalupo et al. 2005, Dijkstra et al. 2006, Verhamme 
et al. 2006, Tasitsiomi 2006). 

In this paper, we will present the implementation and 
validation of Ly-a radiative transfer in LICORICE, a 
dynamics-radiative transfer code, with a special empha- 
sis on the treatment large Hubble flows, which is specific 
to the EoR applications. This work is part of the SKADS^ 
effort (DS2-T1 task) to produce simulated emission maps 
that can be used to optimize the design of the future SKA 
telescope. 

In section 2, we present which physics of the Ly-a 
transfer is included in LICORICE, and justify why some 
aspects of the physics do not need to be included. In sec- 
tion 3 we detail some aspects of the algorithms and accel- 
eration schemes implemented in the code. Three validation 
tests are presented in section 4, comparing the outputs of 
LICORICE against analytic solutions or standard numeri- 
cal results. Finally, in section 5, we apply the code to some 
typical EoR situation to investigate possible fluctuations 
in the Wouthuysen-Field effect. 

2. The code: LICORICE 

LICORICE consists in three parts. A TreeSPH code with 
multiphase modeling of the gas computes the dynamics of 
structure formation (Semelin & Combes, 2002). A contin- 
uum radiative transfer part is added to compute reioniza- 
tion. This part uses a Monte-Carlo approach and is similar 
to CRASH (Maselli et al. 2003). It has the advantage over 
CRASH of using an adaptative grid. LICORICE is cur- 
rently participating in the second part of the Cosmological 
Radiative Transfer Comparison Project (see Iliev et al. 
2006b for the first part of the project). The third part 
is the Ly-a radiative transfer which is described in this 
paper. 



http:/ /www. skads-eu.org 



Semelin, Combes & Baek: Lyman-alpha radiative transfer during the EoR 



3 



2.1. Physics of the Ly-a radiative transfer 

LICORICE implements a Monte-Carlo approach to radia- 
tive transfer: it propagates photons on an adaptative grid. 
Consequently we will describe the physics of Ly-a radia- 
tive transfer from the point of view of a photon traveling 
through the simulation box. 



2.1.1. Computing the optical depth 

A photon propagating through neutral hydrogen has a 
probability P(r) = e~ T of not being scattered after trav- 
eling through an optical depth r from its emission point. 
We consider Ly-a scattering only, so the optical depth can 
be computed as: 



dsdu\\ tihi p(u\\) <J 1/(1 



(1) 



where v is the photon frequency in the global rest frame 
and njji is the local number density of neutral hydrogen, 
it II is the scattering atom velocity along the incoming pho- 
ton's direction in the moving reference frame of the fluid 
, p(u\\) is the normalized probability distribution for mm , 
^macro j g g as macroscopic velocity along the incoming 
photon's direction in the global rest frame, and c is the 
speed of light. Finally <j(i/') is the Ly-a scattering cross 
section of a photon with frequency v' in the atom rest 
frame. 

The function p(u\\) usually results from the thermal 
distribution of the atoms velocity: 



p(u\\) 
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(2) 



Some astrophysical systems have significant velocity gra- 
diants on scales below the best possible resolution of the 
simulations. In these cases it may be relevant to model 
the small scale velocity contribution by adding a turbu- 
lent componant to the thermal velocity dispersion. We did 
not include any turbulent contribution in this paper. 

The exact expression of the Ly-a scattering cross sec- 
tion is given in Peebles (1993). It is well approximated by 
the usual Lorentzian profile: 



l A f ne 
o(v) = fl2 



Ai/ L /2n 



m e c (v - i/ ) 2 + (Av L /2f 



(3) 



where /12 = 0.4162 is the Ly-a oscillator strength, uq = 
2.466 x 10 15 Hz is the line center frequency and Ai/^ = 
9.936 x 10 7 Hz is the natural line width. 

We introduce the dimensionless parameters x, the rel- 
ative frequency shift, and &, the natural to Doppler line 
width ratio: 



v - vq . v th 

x = — with Ai/p — i/q , 

Ai/d c 



Av L 
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(4) 



(5) 



Using these notations, we can write the optical depth 
increment in the gas local rest frame (wf| lacro = 0) as: 



(It = els nni a — ("j x ) ■• 



m e cAvu 



or in cgs units, 



dr = ds n m 6.87 x 1(T 14 ^-j H(b, x) 
where H is the Voigt function defined as: 



T 



H(b,x) = - 



-HOO 



n J -00 - V) 2 + b2 



dy 



(6) 



(7) 



(8) 



To compute the Voigt function, we either use the an- 
alytic fit given by Tasitsiomi (|2006[l . or, in cosmological 
situations where the Hubble flow on a scale of the order 
of the simulation spatial resolution produces a frequency 
shift larger than the Ly-a line width (see section 2.1.3), 
we use the simple approximation: 

H(x, b) ~ max ( , -^-A ■ (9) 



2.1.2. Effect of the n=2 state splitting, dust and 
deuterium. 

A Ly-a photon can excite an hydrogen atom from the 
IS1/2 ground state to either the 2P 1 / 2 or the 2P 3 / 2 state. 
As discussed by Tasitsiomi (2006), the splitting between 
these two n = 2 states is small: only 10 GHz, or 1 km.s -1 . 
If the thermal velocity dispersion is much larger than ~ 1 
km.s -1 (i.c T > 100K), the level splitting is washed out 
from the radiation spectrum after just one scattering. In 
the case of the Ly-a background radiation during the EoR, 
the gas temperature drops to ~ 30K. However the optical 
depth as the photon redshifts from one side to the other 
of the Ly-a line is very high: ~ 8 x 10 5 for the average 
gas density of the universe (fib = 0.045) at z = 9.5 in a 
standard cosmological model. Thus the photon will scatter 
many times off thermal atoms and still the splitting will be 
washed out. Consequently, we did not distinguish between 
the two 2P states in this paper. 

Another issue is the possible reshuffling from 2P to 
2S through collisions with free protons or electrons. An 
atom cannot be excited directly to the 2S state because 
of the dipole selection rule. But if the 2P — > 2S tran- 
sition is induced by a collision, then the atom will de- 
excite through the emission of 2 continuum photons: the 
Ly-a photon is lost. Tasitsiomi (2006, see eq. 26) com- 
putes the probability that an atom in the 2P state will 
reshuffle to 25* before it de-excites normally by emitting 
a Ly-a photon, p ~ 8.5 x 10 _13 n p ( j^r) -0 ' 17 , where n p 
is the proton number density. What can we expect during 
the EoR ? First, we are interested in the Ly-a background 
in the cold, neutral region of the universe, thus n p <C 1. 
Moreover, even if we assume a non negligible ionization 
fraction, p~5x 10~ 13 at z ~ 10 for T — bOK and for 
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the critical density of the universe with Of, = 0.045. This 
probability must be compared with the average number 
of scatterings a photon undergoes as it redshifts through 
the Ly-a line, which is of the order of the optical depth: 
~ 10 6 . We see that only a fraction of 5 x 10~ 7 of the Ly-a 
photons will be degraded into 2 continuum photons before 
they redshift out of the line. We conclude that this process 
is not relevant during the EoR. 



Dust is usually a factor in Ly-a transfer simulations 
since it absorbs Ly-a photons. Hansen and Oh (|2006|) 
study the effect of dust absorption in a multiphase medium 
and show that photons can escape such a medium, while 
they would be absorbed in a more homogeneous single- 
phase medium, for the same total HI column density. In 
the cosmological context we are working in (z > 6), there 
are no observations to help us constrain the dust abun- 
dance and distribution. We can assume that during the 
EoR, dust may be found only around the sources and that 
the IGM is completely dust free. Under this assumption, 
the effect of dust on the Ly-a flux can be modeled with a 
simple escape fraction coefficient. Furthermore, if this co- 
efficient is assumed to be independent of the source, dust 
will not have any effect on the Ly-a flux fluctuations which 
are the focus of this paper. Therefore, we did not include 
the effect of dust. 



Dijkstra et al. p006j) show that the presence of deu- 
terium with an abundance [D/H]= 3 x 10~ 5 leaves a clear 
imprint on the spectrum emerging from a uniform sphere 
of gas with a central source and a total optical depth 
t = 7.3 x 10 5 . Is deuterium relevant to the Ly-a flux 
during the EoR? The deuterium line center is 82 km.s -1 
blueward of the hydrogen line. This is equivalent to the 
redshift of a photon traveling ~ 0.5 comoving Mpc at 
z=9.5. So the first answer is that deuterium may have an 
effect on the Ly-a flux fluctuations only at small scales ( 
< 1 Mpc comoving). However, let us notice that the to- 
tal optical depth for deuterium through the Ly-a line is 
t ~ 20 for an abundance [D/H]= 2 x 10 -5 , and the op- 
tical depth in the wing of the Hydrogen Ly-a line, as the 
photon redshifts from far into the blue to the center of 
the deuterium line happens by chance to be also r ~ 20. 
This means that, while the photon will indeed scatter a 
few times in the deuterium line, it will also have been 
scattered in the hydrogen line wing several times before 
it reaches the frequency range where deuterium scattering 
dominates (a few 10 km.s -1 around the line center). We 
conclude that the presence of deuterium is unlikely to af- 
fect the Ly-a flux fluctuations noticeably. This does not 
mean that the effect of deuterium on Ly-a radiation can- 
not be observed during the EoR. Indeed, in the case of an 
ionizing bubble around a source with a sharp ionization 
front, the continuum spectrum of the source will show a 
Gunn-Peterson trough with a step at the bubble redshift. 
Deuterium should create a secondary small step. 



2.1.3. Dealing with large Hubble flows 

In z ~ 10 cosmological simulations, the Ly-a thermal line 
width is equivalent to the Hubble flow redshift over only a 
few 10 comoving kpc. This scale is usually (much) smaller 
than the size of cells in simulations. So we must be care- 
ful when we compute the optical depth: using a single co- 
moving frequency for the photons throughout a cell would 
result in photons flowing through the line core without 
feeling it. We must actually compute an integral along 
the path of the photon, with a redshifting comoving fre- 
quency, to obtain the correct optical depth. If we consider 
that the expansion velocity between any two points of the 
same cell is non-relativistic (— ~ 0.01 for 20 comoving 
Mpc at z ~ 10), the computation gets easier. Let Ui n be 
the comoving frequency and Xi n the local rest frame value 
of x when the photon enters the cell. 



^™(1 



(10) 



ro is the macroscopic velocity of the gas (uni- 
form inside the cell), and k is the direction of the photon. 
Then, at a given point inside the cell defined by the vector 
rk from the entering point of the photon, the comoving 
frequency is: 



a{v m ) 
v = —r^rVi, 
a(y) 



1 



Hr 



Hr 

f» l ) , (11) 

c 



where H is the Hubble constant at the simulation red- 
shift, and a is the expansion factor. We neglect any vari- 
ation of H during the photon travel and we consider 
non-relativistic expansion velocities. For non-relativistic 
macroscopic velocities of the gas, the corresponding value 
of x writes: 

Hr Vi n 

(12) 



c Aujj 



So, we linearized the relation between r, the current 
path length inside the cell, and the current local rest frame 
value of the x variable. In this approximation, noting x ou t 
the value of x when the photon exits the cell, computing 
the optical depth reduces to computing the integral Voigt 
function Hi n t(x out ), with the following definition for the 
Hi n t function: 



-ffint(x) 



H(x , b)dx 



(13) 



However, in the Monte-Carlo method, to find the lo- 
cation of a scattering event, we need to solve the equation 
Hi n t(x) = A, where A is a constant. This is simplified if 
we can provide an analytic expression for Hi nt (x), with an 
explicit inverse function. It is the case if we use the simple 
approximation of H(x, a) given in eq. [9] It involves the 
erf(x) function for which we are using an approximation 
which has an explicit inverse function. 

Modeling the effects of expansion only by a redshift 
computed from a radial dilatation is the usual approxi- 
mation for Ly-a radiative transfer codes. Other effects of 
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the variation of the expansion factor a, such as variation 
of the average density during the photon flight time, are 
usualy ignored. However, let us emphasize that adding the 
expansion velocity to the other types of velocities (macro 
or microscopic) to compute Doppler shifts, either dur- 
ing scattering events or to compute local rest frame val- 
ues of x, works only if all velocities are non-relativistic. 
However, we want to study the Ly-a flux during the 
EoR, and at z ~ 10, a photon emitted just below the 
Ly-/3 frequency will travel <~ 350 comoving Mpc and be 
redshifted to Ly-a by a Hubble flow velocity such that 
vr/c ~ — ~ — ~ 0.15. In this case, the second order er- 
rors in computing the redshift are not completely negligi- 
ble. But, what is more important, neglecting the variations 
of a along the photons path, in computing gas densities for 
example, produces first order error in — . Consequently, 
we should limit ourselves to ~ 100 comoving Mpc boxes. 
In this work we do neglect variations of a, except for the 
cosmological redshift, and in most cases we use simulation 
boxes smaller than ~ 30 comoving Mpc. The full effect of 
expansion will be introduced in the future to handle larger 
boxes. 



2.2. Scattering off hydrogen atoms 

In an expanding universe the natural variable is the co- 
moving frame frequency of the photon. However, as long 
as the Hubble flow velocities are non-relativistic, we can 
use the frequency in the rest frame of the zero-coordinate 
point of the simulation, hereafter named global frame, and 
treat the expansion as a simple radial dilatation with a 
vjj = Hr velocity field to be added to the peculiar veloc- 
ities. 

In the rest frame of the atom, we will consider the 
scattering to be resonant. The effect of the recoil, which 
would change the frequency of the photon by transferring 
part of its momentum to the atom, has been shown to 
be negligible in astrophysical situations by several authors 
(Zheng & Miralda-Escude l2"0"02"l or Tasitsiomi lStM]) . In the 
global frame, however, due to the various contributions 
to the atoms velocity, the frequency v of the photon will 
change. 

There are three main contributions to the atom veloc- 
ity: the thermal motion, the macroscopic peculiar motion 
and the Hubble flow. When the photon scatters off an 
atom, we first compute the frequency in the atom rest 
frame: 



^atom 



«/(!-( 



ki) , 



(14) 



where u is the thermal velocity of the atom, and is 
the direction of the incoming photon. Let us split u into 
U||, the component parallel to the incoming photon direc- 
tion, and ui, the perpendicular component. obeys the 
distribution: 



while each of the 2 components of obeys: 



P*(y) = 



'KVth 



(16) 
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nH{b, x) [x - y) 2 + b 2 



(15) 



Then we compute the direction of the photon after scat- 
tering. As described by Tasitsiomi (|2006|) the scattering 
phase function depends on the excitation state and on 
whether the photon scatters in the wing or in the core of 
the line. Several authors (see e.g. Cantalupo et al. 2005]or 
Verhamme et al. 120061 ) have shown that, for high opti- 
cal depth media, the shape of the phase function does not 
alter the results of the simulation. In this work, we use 
isotropic scattering. Using the new direction of the pho- 
ton, we then recompute the frequency in the global frame 
using a non-relativistic Doppler effect. 



3. Ly-a radiative transfer: numerical methods 

3.1. Gas density and velocity field 

LICORICE is meant to use results from dynamical simu- 
lations done using the TreeSPH part of the code (Semelin 
and Combes 2002). It uses the same Tree structure as the 
dynamical code to build an adaptative grid. The grid is 
build in such a way that each cell contains between 1 and 
iV max gas particles. Values from 1 to 30 are commonly 
used for N max . The density and velocity fields are then 
interpolated from the particle distribution. The cells limit 
the resolution of the simulation only by having uniform 
dynamical properties. Ly-a transfer inside one cell is still 
computed exactly under the assumption that the expan- 
sion velocity on the scale of the cell is non-relativistic. 

3.2. Monte-Carlo method 

Using the Monte-Carlo approach, we send individual pho- 
ton from the source and follow them from scattering to 
scattering and from grid cell to grid cell, until they exit 
the simulation box. After one event (emission or scatter- 
ing), the algorithm is the following: 

— Step 1: Compute the photon global frame frequency, 
either from the scattering atom rest frame frequency, 
or from the source spectrum. 

— Step 2: Draw the new photon direction (isotropically 
in this work). 

— Step 3: Draw a variable p from a uniform distribution 
between [0, 1]. The photon will travel an optical depth 
r = — ln(p) to the next scattering event. 

— Step 4: Increment optical depth with current cell con- 
tribution. Determine if scattering occurs in this cell, if 
yes go to step 5, if no, pass on to next cell and repeat 
step 4. 

— Step 5: Draw scattering atom thermal velocity, and 
compute frequency in the scattering atom rest frame. 
Go back to step 1. 
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3.3. Acceleration scheme 



To generate random variable following eq . 1 1 51 distribution, 
we use the method by Zheng & Miralda-Escude (|2002p . 
These authors introduce a parameter uq for which an op- 
timal value is needed. We use the following empirical fit: 

M = 1.85-log(6)/6.73 + ln(ln(x)) for x > 3 , (17) 

and Mo = otherwise. This has been determined from 
a systematic numerical optimization, and works well for 
10~ 4 < b < 10 _1 . For large values of x, quite common 
in cosmological simulations, drawing a value from Pi is 
still slow, even with Zheng & Miralda-Escude method. In 
this case, however, scattering will most likely occur in the 
wing of the line. Consequently, for x > 10, we revert to 
uq = but we truncate the distribution to the limited 
range [—3, 3]. 

We use the core-skipping acceleration scheme (see 
Avery & House 1968 or Ahn, Lee and Lee 2002 for first 
applications). In this scheme, we choose a core value x c for 
the variable x. In the regime x < x c the medium must be 
thick. If so, the photon will scatter many times over a very 
short distance. Only when x gets larger than x c (scatter- 
ing in the wing by a fast moving atom) will the medium 
become transparent and will the photon travel a long dis- 
tance. The idea is to ignore the insignificant core scatter- 
ings: every time the photon enters the core (x < x c ), it 
leaves again immediately by scattering off an atom with 
a thermal velocity u at > x c Vth- We use the detailed pre- 
scriptions given by Tasitsiomi (|2006[) on choosing Xq as a 
function of &to, where b is defined in eq. 5 and tq is the op- 
tical depth at the line center . The core-skipping scheme 
works well when the expected output of the code is an 
emerging spectrum. We will show that it also works for 
computing the fluctuations of the local scattering rate. 

4. Validation tests 

LICORICE being a complex, multipurpose code, we take 
care of validating each part separately. We present here 
validation tests for the Ly-a part, against analytical solu- 
tions or standard numerical setups. 

4.1. Static homogeneous sphere: emerging spectrum 

A classical validation test for Ly-a codes is the emerg- 
ing spectrum for a monochromatic source in the middle 
plane of a static homogeneous slab of gas. The main rea- 
son is that Neufeld (|1990|) gives an analytic solution for 
the emerging spectrum in the case of an extremely thick 
system. However Dijkstra et al. f)2006[) provide a new an- 
alytic expression in the case of a uniform spherical cloud 
of gas: 
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Fig. 1. Emerging spectrum for a uniform and spherical 
cloud of gas at T= 10 K and several values of To, the 
total optical depth at line center. Numerical histograms 
are computed with 10 5 photons. The analytic solution is 
given in Dijkstra et al. (|2006|> . 



Since this geometry is more relevant in a cosmological con- 
text, we present results for this case (we also performed 
the Neufeld test and also found a good agreement). Fig. 
[T] shows the comparison between numerical and analytic 
spectra, for a static spherical cloud of gas with temper- 
ature T = 10K and optical depths at line center from 
center to edge equal to 10 5 , 10 6 and 10 7 . For each run, 
10 5 photons are injected at the center of the cloud with 
frequency vq. The emerging spectrum J(x) is obtained 
by computing a normalized histogram of the frequencies 
of the photons as they leave the cloud. The agreement is 
good in all cases, although it gets better and better as To 
increases. Indeed in the case r = 10 5 , we have aro = 1500 
which is close to the lower limit for an extremely thick 
medium. 



4.2. Expanding homogeneous sphere: emerging 
spectrum 

The second test is also becoming a classic. It was first 
performed by Zheng and Miralda-Escude (2002), followed 
by several authors. Verhamme et al. (2006]) present espe- 
cially detailed results. The system is a uniform sphere of 
gas expanding (or contracting) with an Hubble-like ve- 
locity field: the radial velocity is proportional to the ra- 
dius. We use the same physical condition as Zheng and 
Miralda-Escude (|2002|) : a temperature of 20 000K and a 
radial velocity of 200 km.s -1 at the edge of the cloud. We 
use three different column densities from the center to the 
edge of the cloud: 7V H equals 2. x 10 18 cm~ 2 , 2. x 10 19 
cm" 2 or 2. x 10 20 cm" 2 (that is r = 8.3 x 10 4 , 8.3 x 10 5 , 
or 8.3 x 10 7 ). We do not run the tests for a contracting 
cloud, and we consider only a central point source emitting 
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Fig. 2. Emerging spectrum for a spherical, uniform and 
expanding cloud of gas at T= 20000 K and several values 
of Nhi, the column density from the center to the edge 
of the cloud. The emission is at Ly-a frequency from the 
central point of the cloud. The radial velocity of the gas 
is proportional to the radius and equals 200 km.s -1 at 
the edge of the cloud. The peak blueward of the Ly-a 
frequency is completely suppressed by the expansion. 



Fig. 3. Mean intensity spectra at various radii for a Ly-a 
monochromatic source in a uniform expanding medium at 
T = OK. The numerical result is compared to the ana- 
lytic solution given, in the diffusion regime, by Loeb and 
Rybicki (1999). 



at Ly-a frequency, not the case of uniform emissivity. An 
expanding cloud with a central point source is the most 
relevant to the EoR. The emerging spectra are shown in 
fig. [H We find results very similar to Zheng and Miralda- 
Escude ((20021) or Verhamme et al. ([2006]) . The peak blue- 
ward of the Ly-a frequency is completely suppressed by 
the expansion (not shown in fig. [2]). For the case of a con- 
tracting cloud, the red peak would be suppressed 

4.3. Expanding homogeneous sphere at T=0K: mean 
intensity field 

Our last validation test relies on an analytic solution 
given by Loeb and Rybicki (1999) . The setup is again 
an expanding uniform hydrogen medium with a central 
monochromatic source at Ly-a frequency. The velocity 
field is a Hubble flow again, but this time, the gas tem- 
perature is T = OK. Loeb and Rybicki introduce the di- 
mensionless variable v = 1/0 ~ v , where is the comoving 
frequency shift from Ly-a at which the optical depth to 
infinity equals 1, and f = -7-, where r* is the proper ra- 
dius at which the Doppler shift from the source due to 
the Hubble expansion equals v*. Loeb and Rybicki give 
an analytic expression for the corresponding dimension- 
less mean intensity J, valid in the diffusion regime: 
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Loeb and Rybicki l|1999[) show the comparison between 
the analytic solution and the numerical solution given by a 



dedicated Monte-Carlo code. Tasitsiomi (|2006p runs this 
test with his general purpose code. Our results are pre- 
sented in fig. [3] and [4] They are very similar to those 
of Loeb and Rybicki and Tasitsiomi: the numerical re- 
sults are close to the analytic solution where the diffusion 
regime is valid. However, where photons enter the free 
streaming regime (fo = 1 in fig. [3] and log^o = 0.5 in 
fig. [4j, the numerical solution diverge from the analytic 
solution which becomes invalid. 

Here are some details on how we ran this test. We used 
a temperature T = 2K for computing the optical depth, 
which is singular for T = OK, but we used a zero thermal 
velocity for the atoms in all scattering events. We used 
a simulation box holding a sphere of gas of radius lOr*. 
Since the temperature is not zero in all respects, it rein- 
troduces a dimension in the problem. Very close to the 
source the thermal speed is larger than the expansion ve- 
locity and the numerical behaviour should diverge from 
the analytic solution. For reference, we used a proper ex- 
pansion velocity of 200 km.s -1 at r*. In our setup, the 
thermal and expansion velocity are of the same order for 
r ~ lO -3 ^. 

The radii range covered by this test is quite large, from 
10 -3 rv to lOr*. Using 10 12 cells for the radiative transfer 
grid is not really an option. Actually we used only 16 3 
cells, but we took advantage of the integral scheme de- 
scribed in equation 12 and 13 for computing exactly the 
optical depth between two points in the same cell. This 
approach is obviously validated by the good agreement 
shown in fig [4] at small radii. 
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Fig. 4. Mean intensity profile at various frequencies for 
a Ly-a monochromatic source in a uniform expanding 
medium at T — OK. The numerical result is compared 
to the analytic solution given, in the diffusion regime, by 
Loeb and Rybicki (1999). 

5. Ly-a during the EoR 

5.1. The role of Ly-a photons in the 21 cm emission 

The 21-cm signal can be seen in emission or absorption 
against the CMB. The differential brightness temperature 
observed at redshift z=0 is: 

ST b = Ts ~^™ B {l-e- T21 ), (20) 

where Tg is the neutral hydrogen spin temperature, Tqmb 
is the CMB radiation blackbody temperature at redshift 
z, and T21 is the 21-cm line optical depth. The value of t 2 i 
is given, among others, by Tozzi et al. (2000) or Furlanetto 
et al. (2006a). Injecting the formula for T21 -C 1, the dif- 
ferential brightness temperature can be written: 

5T b ~ 9. xm (1 + 5)(1 + z)i Ts ~J CMB mK , (21) 

where 5 is the local overdensity at redshift z, and xm is 
the neutral hydrogen fraction. This value of ST b is for 
a flat ACDM model with h = 0.7, and Q, b = 0.044. It 
changes by ± 0.5 mK when Cl m varies from 0.25 to 0.3. 
The value of the spin temperature Ts is the result of three 
competing processes. The absorption/reemission of CMB 
photons tends to bring Ts to X"cmb over a time scale under 
10 5 years during the EoR (Tozzi et al., 2000). As we have 
seen in the Introduction section, both collisions between 
hydrogen atoms and the pumping by Ly-a photons, also 
known as the Wouthuysen-Field effect (Woutuysen 1952, 
Field 1958), tends, instead, to couple Ts to the kinetic 
temperature of the gas. As a result, the spin temperature 
can be written (Furlanetto ct al. 2006a): 



Tg 1 ^ T cmb + x cT K 1 + x a T c 1 w . th Tc ^ Tr ^ 
1 x c -\- Xq, 

where x c and x q, are the coupling coefficients respectively 
for collisions and Ly-a pumping, and Tq is the effective 
color temperature of the UV radiation field (see Furlanetto 
et al. 2006a) . The coefficient x a , which is the focus of this 
work, can be explicitly written as: 

27Ai Tcmb 

where T* = 0.068K, A w = 2.85 x lO" 1 ^- 1 is the sponta- 
neous emission factor of the 21-cm transition, and P a is 
the number of scatterings of Ly-a photons per atom per 
second. Now usually come two approximations that we 
will not make in this paper. The first is that P a is con- 
sidered proportional to J{y a ), the angle averaged specific 
intensity at the local Ly-a frequency, neglecting the contri- 
bution of wing absorptions. We will see that this approx- 
imation is valid provided J{v) itself has been computed 
taking into account wing absorptions. The second, more 
drastic approximation, is to evaluate J{v) without per- 
forming the full radiative transfer computation. Actually, 
to our knowledge, all numerical simulations of 21-cm emis- 
sion consider a uniform value of J{v a ). However, Barkana 
and Loeb (2005) have shown that several factors induce 
fluctuations in J{v a ): the 1/r 2 scaling of the flux which 
magnifies the Poisson noise in the source distribution, the 
clustering of the sources, and the contribution of higher 
Ly-a series photons (also studied in detail by Pritchard 
and Furlanetto 2006). They predict the power spectrum 
of the 21-cm brightness temperature due to the fluctua- 
tion in J(v a ). Although a vast improvement over using a 
uniform Ly-a flux, they are still neglecting radiative trans- 
fer effects: they assume that photons arc freely streaming 
until they redshift to the local Ly-a frequency. We will 
show that this assumption breaks down at scales smaller 
that ~ 10 comoving Mpc. 

Potentially, pumping the upper excitation level is not 
the only way in which Ly-a photons can influence the 
21-cm emission: they also heat up the gas. This heating 
mechanism was first thought to be efficient by Madau et al. 
(1997). However, Chen & Miralda-Escude (2004), taking 
into account the atoms thermal velocity distribution which 
had been neglected by Madau et al., found a much smaller, 
actually negligible heating rate. Furlanetto & Pritchard 
(2006c), taking into account the effect of higher Lyman 
series photons confirmed this result and found, for typical 
EoR conditions, the Ly-a heating rate to be 140 smaller 
than the heating rate from X-rays. Chuzhoy & Shapiro 
(2007a) recently challenged this result, adding in particu- 
lar the effect of the deuterium Ly — (3 resonance line, but 
the strength of this effect is not yet completely probed. 
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Fig. 5. Profiles of the scattering rate per atom,i- > Q ,, in the 
Ly-a line, around a central source with a continuous flat 
spectrum in an homogeneous medium of neutral hydrogen 
at 30K and the average baryon density of the universe at 
z ~ 10. 

5.2. P a profiles for spherically symmetric 
configurations 

In the next three cases, we consider a central source in 
a spherically symmetric medium of neutral hydrogen at 
T K = 30K, and z ~ 10, which is typical of the EoR . The 
main difference with the setup of the tests in section 4.3 is 
that the source emits a continuous flat spectrum. We deal 
only with photons between Ly-a and Ly-/3 frequencies. 

5.2.1. Homogeneous medium 

First, we consider the case where the gas is homogeneous 
with a density equal to the average density of the universe 
at z ~ 10. We ran the simulation in a 600 comoving Mpc 
box with the source in the center of the box. We present 
the results up to a radius of 100 comoving Mpc: at larger 
radii the effect of the box boundaries (no retrodiffusion) 
and the variation of a during the photon flight (see section 
2.1.3) alter the validity of the simulation. Fig. [5] presents 
the radial profile of the scattering rate per atom, P a , in 
different cases. Using the real Voigt line profile, we can 
see a deviation at small scale from the simple 1 /r 2 profile 
expected if photons stream freely from the source until 
they redshift to the local Ly-a frequency. Here is why. Let 
us consider a photon emitted above the Ly-a frequency. 
As it travels away from the source, it is redshifted toward 
the local Ly-a frequency. Because of the contribution of 
the wings of the Voigt profile to the optical depth, it has 
a probability to scatter before reaching the local Ly-a fre- 
quency On average, it will scatter for the first time, when 
the optical depth along its path reaches 1. This is achieved 
when the photon is redshifted to the frequency ~ v a + v* 
(see section 4.3), which occurs at a distance r* before the 
location where it would reach the local Ly-a frequency 



(quantities defined in Loeb & Rybicki 1999). In our setup, 
typical of the EoR, r+ ~ 10 comoving Mpc. After the 
photons scatter for the first time they change direction. 
Consequently, the location where a photon actually en- 
ters the core of the local Ly-a line has a probability to be 
anywhere within ~ 10 comoving Mpc of the location deter- 
mined by free streaming alone. This is not crucial at large 
scales where the expected 1/r 2 profile is recovered, but it 
creates a steeper profile at small scales (exponent ~ —2.3 
in our setup). To validate our interpretation, we computed 
the transfer with a modified line profile: the core of the 
line is unchanged, but the wings are set to zero. As can 
be seen in Fig. [51 we then recover the 1/r 2 profile. After 
this paper was submitted, Chuzhoy and Zheng (2007b) 
posted a paper with a similar result. They considered a 
very similar setup and computed the transfer with a sim- 
ple Monte Carlo code which is 1-D (spherical symmetry) 
and uses a simplified line profile but does include Ly-a 
photons locally injected by cascades from upper Lyman 
series lines. They find the same steepening of the scatter- 
ing rate profile at short scales. However, they show that 
photons injected from upper Lyman series lines are much 
less sensitive to radiative transfer effect: they follow the 
1/r 2 profile more closely. Consequently the discrepancy 
between the full radiative transfer computation and the 
simple 1 /r 2 evaluation is somewhat reduced. 

Finally we checked the effect of using the core-skipping 
acceleration scheme for evaluating the fluctuations of P a : 
since all core absorptions are avoided, it could have modi- 
fied the spatial P a fluctuation map. As can be seen in Fig. 
[5j it is not the case. Although the wings modify the loca- 
tion where photons enter the core, the scattering number 
is still dominated by the close-to-the-core contribution. 

5.2.2. Central clump 

Until now, the Ly-a scattering rate has been evaluated 
during the EoR in an homogeneous medium only (Barkana 
and Loeb 2006, Pritchard and Furlanetto 2006). With a 
general purpose 3D code such as LICORICE, we can inves- 
tigate the impact of density fluctuations in the gas on the 
Ly-a scattering rate per baryon. First we choose a very 
simple setup: we consider a box of size 64 comoving Mpc, 
with a central homogeneous spherical clump of gas of ra- 
dius 1 comoving Mpc. The clump is 64 times denser than 
the surrounding medium, which has the average baryon 
density of the universe at the redshift of the simulation 
(z ~ 10). The source is in the center of the clump and has 
a continuous flat spectrum. The radial profile of the Ly-a 
scattering rate per atom, P a , is shown in Fig. [5] The main 
feature is a depletion of P a in the low density medium just 
outside the clump. Indeed, photons that should redshift 
to the Ly-a frequency in this region, have first to travel 
through the high density region where they have an en- 
hanced probability to be scattered in the wing of an atom, 
and redirected to redshift to local Ly-a while still inside 
the core. In other words, the enhanced wing scatterings 
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Fig. 6. Profile of the scattering rate per atom, P a , in the 
Ly-a line, around a central source with a continuous flat 
spectrum, inside an overdense spherical clump of gas (64 
times the density of the surrounding medium) of radius 1 
comoving Mpc, 



in the clump draw core scatterings from surrounding re- 
gions to the clump itself. This is somewhat similar to the 
usual shadowing effect in radiative transfer, although in 
this case the process occur in the frequency space so the 
shadowing can be seen even in a spherical configuration. 
Obviously, in addition to the l/r 2 decline of the flux, and 
the fluctuations in the source distribution, we can expect a 
new source of fluctuations for P a : the density fluctuations 
of the intergalactic medium. 

5.2.3. Isothermal density profile 

Instead of the sharp density transition between the clump 
and the surrounding medium we now consider an isother- 
mal density profile (p ~ l/^ 2 ) around the source. The 
isothermal profile connects to the surrounding homoge- 
neous medium at a radius of 10 comoving Mpc. All other 
parameters of the simulation are the same as in the pre- 
vious setup. The radial profile of the Ly-a scattering rate 
per atom, P a , is shown in Fig. [7j Since the average density 
is higher than in the other setups, photons scatter many 
times and even with a acceleration method, we used only 
2 x 10 4 photons. To avoid a very high noise level at large 
radii we used an adaptative resolution: high in the cen- 
ter, lower in the outer regions. We can see on the figure 
that, inside the region with an isothermal density pro- 
file, P a closely matches a l/r 3 profile. It reverts to l/r 2 
in the homogeneous medium outside the 10 Mpc radius. 
This shows that the brightness temperature fluctuations 
of the 21-cm emission may be stronger than previously 
estimated at small scales, at least during the early EoR 
when P a fluctuations are meaningful. 

Is this setup more relevant to the prediction of P a fluc- 
tuations during the EoR than a uniform medium ? Yes in 




Radius (comoving Mpc) 

Fig. 7. Profile of the scattering rate per atom P a in the 
Ly-a line, around a central source with a continuous flat 
spectrum. Inside a sphere of radius 10 comoving Mpc, the 
gas density field is an isothermal sphere, outside it is ho- 
mogeneous at the average baryon density of the universe 
(z ~ 10). The two regions connect smoothly. A l/r 3 profile 
is plotted for comparison. 



the sense that the medium should obviously be denser 
closer to the source. However the spherical symmetry and 
the specific profile used here are oversimplified: the actual 
intergalactic medium in the EoR is not in an equilibrium 
configuration, especially during this early period of sources 
formation, and the spherical symmetry is broken by the 
presence of filaments. In the next section we investigate 
how filaments modify P a . 

5.3. P a map for an axisymmetric configuration 

We consider now an axisymmetric density field for the gas: 
inside a cylinder of radius 1 comoving Mpc the density is 
64 times the critical baryon density, outside it is equal to 
the critical baryon density. The source is located on the 
symmetry axis. All other parameters are identical to the 
previous setup. This setup is a simplified representation 
of the real case of source formation during the EoR at the 
intersection of several filaments, where the filaments have 
a density profile. The P a contour map is shown in Fig. 
[H The shaded area represents the filament, and the con- 
tours are equally spaced on a logarithmic scale (two per 
decade). The map is integrated over a variation of ir of 
the angular variable. About 5.10 6 photons where used for 
this simulation. There is some boundary effect due to the 
finite size of the simulation box: photons that would scat- 
ter just outside the simulation box and possibly reenter 
the box are lost instead. This affects a few Mpc near the 
boundary. We can see a sharp depletion in the number of 
scatterings per atom inside the filament. We have checked 
that a smaller density contrast creates a smaller deple- 
tion, as expected. Once again, we see that the gas density 
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Fig. 8. Contours of the scattering rate per atom P a in 
the Ly-a line, around a central source with a continuous 
flat spectrum. The source is located inside a cylindrical 
homogeneous filament of gas (shaded zone in the figure) 
with a density 64 times the density of the surrounding 
medium which is at the average baryon density of the 
universe (z ~ 10.). The contours are equally spaced on a 
logarithmic scale with a step of VTO. 

field fluctuations induce fluctuations in P a . However, the 
weaker Wouthuysen-Field effect in denser regions may be 
balanced by the greater coupling due to collisions (which is 
proportional to the density) . What may be more relevant 
to the future observations than what occurs at small scales 
inside the filaments, is that the presence of the filaments 
modifies the shape of the contours in the low density sur- 
rounding medium: they are not spherical but oblate. With 
our density contrast, the axis ratio is about 2 in the 5-10 
Mpc range and decreases at large distances. Fig. [9] shows 
the normalized P a profiles along filament axis and in the 
perpendicular plane containing the source. While inside 
the filament (distance smaller than 1 comoving Mpc), the 
profiles match. At larger distances, fig. [9] quantify the P a 
ratio between the two regions: it reaches a maximum value 
of ~ 100 at a distance of ~ 5Mpc, 

6. Conclusions 

The main goal of this paper was to investigate a source 
of fluctuations in the brightness temperature of the 21-cm 
emission during the EoR usually neglected in numerical 
simulations. The Wouthuysen-Field effect, which couples 
the hydrogen spin temperature to the kinetic tempera- 
ture of the gas, is regulated by the crucial parameter P a : 
the number of scatterings per atom per second. While all 
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Fig. 9. Scattering rate per atom P a in the Ly-a line along 
the axis of a filament and in a perpendicular plane. A cen- 
tral source with a continuous flat spectrum is located at 
the intersection of the filament and the plane. The fila- 
ment is cylindrical with radius 1 comoving Mpc and ho- 
mogeneous with a density 64 times the density of the sur- 
rounding medium which is at the average baryon density 
of the universe (z ~ 10.). 

previous simulations of the 21-cm emission used a uni- 
form value for P a , Loeb and Barkana (2005) have shown, 
in a simple theoretical framework, that several sources of 
fluctuations in P a can modify the power spectrum of the 
21-cm emission. We studied how a full 3D radiative trans- 
fer treatment of the Ly-a line in a cosmological context 
modifies the picture given by Loeb and Barkana. 

The first step was to implement and validate the Ly-a 
radiative transfer in LICORICE. We used a Monte-Carlo 
approach, and implemented an algorithm and acceleration 
schemes similar to those of other existing codes. We dis- 
carded physical processes such as recoil or deuterium con- 
tribution, which have negligible effect for the Ly-a trans- 
fer during the EoR. On the other hand we took care to 
compute the optical depth accurately in an expanding cos- 
mological medium, without any resolution effects. We pre- 
sented three validation tests for the code. The first two are 
classical setups: a monochromatic source in a static uni- 
form sphere of gas, and in an expanding sphere of gas. The 
agreement for the emerging spectrum with analytic solu- 
tions and results by other authors is good. The third test, 
the mean intensity map for a monochromatic source in an 
expanding sphere of gas at T = OK, focuses on a quantity 
more closely related to P a . The comparison with the ana- 
lytic solution provided by Loeb and Rybicki is good where 
the analytic solution is valid: in the diffusion regime. This 
set of tests strongly suggests that the LICORICE is valid. 

Barkana and Loeb (2005) show that several factors 
contribute to the fluctuations in the local Ly-a flux: the 
1/r 2 scaling of the flux, and the Poisson noise and clus- 
tering in the sources distribution. They assumed, however 
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that photons rcdshift freely from the source until they 
reach the local Ly-a frequency, and only then scatter off 
hydrogen atoms. In other words, they neglected wing scat- 
tering. We computed the P a profile for a source with a flat 
continuous spectrum in a uniform expanding medium. We 
showed that the effect of taking into account the wing 
scatterings in a full radiative transfer code is to steepen 
the profile at small scales to a ~ 1/r 2 - 3 profile (with our 
choice of parameters). At large scales, or when wings are 
suppressed, we recover the 1/r 2 profile. Thus we may ex- 
pect, at small scales, more power in the 21-cm emission 
than predicted by Barkana and Loeb. But a yet stronger 
effect was obtained when we introduced fluctuations in the 
density of the gas surrounding the source. We investigated 
a central clump setup and an isothermal profile. In both 
cases we observed alterations in the P a profile. In the case 
of an isothermal density profile we found a ~ 1/r 3 pro- 
file for P a . This also suggests stronger fluctuations than 
predicted by Barkana and Loeb, but, once again, mainly 
at small scales since the as at z — 10 is only ~ 0.1: large 
scale structures had not enough time to grow yet. 

Finally, we tried to create a more realistic situation by 
placing a flat spectrum source inside a filament of over- 
dense gas. We observed a sharp depletion of P a inside the 
filament. The photons scatter out of the filament before 
they reach the core of the line. This fluctuation inside the 
filament is once again a small scale feature, and may be 
difficult to catch with LOFAR or SKA. However, the pres- 
ence of the filament also produced oblate contours for P a 
at larger scales (> 10 comoving Mpc), in the surround- 
ing medium. This effect may be more within reach of the 
resolution of these instruments. 

The increased fluctuations of P a due to radiative trans- 
fer effects and to the inhomogeneous distribution of the 
gas translate linearly into fluctuations of the x a coeffi- 
cient. On one hand, we expect these fluctuations to be 
globally significant and produce brightness temperature 
fluctuations only as long ~ 1, i. e. as long as the 

coupling does not saturate to T s — Tk ■ This occurs in the 
early phase of reionization. It is not possible to be much 
more specific in terms of redshift because the Ly-a pump- 
ing efficiency depends strongly on the model for the source 
type and formation history. On the other hand we showed 
that density fluctuations in the gas can create fluctuations 
of P a of a factor greater than 10. The depleted regions 
will fill up only when the amount of young sources gets 
more than 10 times larger. Cosmological simulations sug- 
gest that this corrresponds to a change of redshift between 
1 and 2 around redshift 10. Consequently we expect that 
the fluctuations of P a due to inhomogeneous gas will lead 
to a longer survival of depleted regions where T$ remains 
coupled to Tcmb- One may argue that the strongest fluc- 
tuations of the gas density are located around the sources 
and are ionized very early. However numerical simulations 
have shown that the ionization front is not spherical: it 
is trapped in the high density regions such as filaments 
pointing to the source, where reionization is much delayed 



(see for example, Gnedin 2000). So, these filaments should 
be able to play their role in creating P a fluctuations. 

We have not investigated the effect on P a of an 
anisotropic peculiar velocity field around the source. In 
principle, it would also induce non spherical contours. 
However, cosmological simulations in a 20 Mpc box sug- 
gest peculiar velocities of the order of 100 km.s -1 during 
the EoR, when the Hubble constant is H(z=10) ~ 1000 
km.s - .Mpc -1 . Moreover, only velocity differences will al- 
ter P a contours, which should be smaller than 100 km.s -1 
at large scales. We estimate that the impact of the velocity 
field of the gas is smaller than the impact of the density 
fluctuations. However, LICORICE fully implements the 
effects of the gas peculiar velocity, and it will be taken 
into account in the future simulations of a cosmological 
box. 

Another process will have to be included in the future: 
the effect of higher Lyman series lines. Barkana & Loeb 
(2005) and Pritchard and Furlanetto (2006) have shown 
that these lines, having horizons closer to the source than 
the Ly-a line, add to the power of the P a fluctuations 
close to the sources. In a forthcoming paper, we will ap- 
ply LICORICE to a cosmological field during the EoR 
using simulation outputs from the HORIZON projectd, 
and compute the resulting 21-cm brightness temperature 
map. If, as we believe, we find significant modification of 
the predictions, we will include higher Lyman series lines. 
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